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Abstract.  In  a  biological  cell,  cellular  functions  and  the  genetic  regula¬ 
tory  apparatus  are  implemented  and  controlled  by  a  network  of  chemical 
reactions  in  which  regulatory  proteins  can  control  genes  that  produce 
other  regulators,  which  in  turn  control  other  genes.  Further,  the  feed¬ 
back  pathways  appear  to  incorporate  switches  that  result  in  changes  in 
the  dynamic  behavior  of  the  cell.  This  paper  describes  a  hybrid  systems 
approach  to  modeling  the  intra-cellular  network  using  continuous  differ¬ 
ential  equations  to  model  the  feedback  mechanisms  and  mode-switching 
to  describe  the  changes  in  the  underlying  dynamics.  We  use  two  case 
studies  to  illustrate  a  modular  approach  to  modeling  such  networks  and 
describe  the  architectural  and  behavioral  hierarchy  in  the  underlying 
models.  We  describe  these  models  using  Charon  [2],  a  language  that 
allows  formal  description  of  hybrid  systems.  We  provide  preliminary  sim¬ 
ulation  results  that  demonstrate  how  our  approach  can  help  biologists 
in  their  analysis  of  noisy  genetic  circuits.  Finally  we  describe  our  agenda 
for  future  work  that  includes  the  development  of  models  and  simulation 
for  stochastic  hybrid  systems.1 


1  Introduction 

In  order  to  survive,  organisms  continuously  monitor  their  surroundings  and,  if 
necessary,  adjust  traffic  through  simple  or  complex  combinations  of  genetic  and 
metabolic  networks  to  respond  to  alterations  in  local  conditions.  Local  condi¬ 
tions  include  both  the  physical  environment,  for  example,  temperature  (the  heat 
and  cold  shock  response),  nutrient  and  energy  source  concentrations  (the  strin¬ 
gent  response),  light  (circadian  rhythms),  cell  density  (quorum  sensing  response) 
as  well  as  the  molecular  environment  of  individual  regulatory  components.  Ex¬ 
amples  of  the  latter  include  intracellular  concentrations  of  transcription  factors 
and  allosteric  effectors.  The  availability  of  complete  genomic  information  for  a 
wide  variety  of  organisms  and  the  consequent  attention  on  proteomics  has  dra¬ 
matically  increased  the  number  of  systems  and  components  of  systems  that  are 
involved  in  these  sensing  and  responding  activities  [4,10].  Understanding  how 
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these  biological  systems  are  integrated  and  regulated  and  how  the  regulation 
may  be  influenced,  possibly  for  therapeutic  purposes,  remains  a  significant  chal¬ 
lenge. 

In  this  paper  we  model  and  simulate  examples  of  genetic  and  metabolic  net¬ 
works  using  a  hybrid  systems  approach  that  combines  concepts  and  tools  from 
control  theory  and  computer  science.  First  we  analyze  a  previously  published 
plasmid-based  genetic  network  that  was  designed  and  synthesized  using  three 
repressor  transcription  factors  where  one  repressor  negatively  regulates  the  pro¬ 
duction  of  a  subsequent  repressor  [7].  Then  we  model  a  biologically  important 
genetic  network  that  controls  the  quorum  sensing  response,  an  adaptive  response 
of  certain  gram  negative  bacteria  to  local  population  density  [13,17] .  The  quorum 
sensing  response  controls  the  luminescent  behavior  in  certain  strains  of  Vibrio 
which  has  been  linked  to  the  normal  development  of  the  bacterial  host  [18]  as 
well  as  to  medically  important  phenomena  such  as  biofilm  formation  by  Pseu¬ 
domonas  aerugenosa,  an  organism  that  can  cause  overwhelming  pneumonia  and 
septic  shock  [11,20]. 

2  Modeling 

The  genetic  circuits  and  biomolecular  networks  considered  here  and  elsewhere 
are  remarkably  similar  to  hybrid  systems  encountered  in  engineering,  for  exam¬ 
ple  embedded  systems.  In  particular,  it  is  worth  noting  the  following  three  key 
features: 

Concurrency  and  communication.  At  the  intra-cellular  level,  proteins  and 
mRNAs  are  agents  communicating  with  each  other  and  influencing  each 
other’s  behavior.  At  the  inter-cellular  level,  cells  can  be  viewed  as  networked 
agents  interacting  with  each  other  via  different  communication  mechanisms. 
Discrete  and  continuous  behaviors.  At  the  lowest  level,  the  evolution  of  en¬ 
tities  such  as  proteins  can  be  described  by  differential  equations.  Discreteness 
arises  in  two  ways.  First,  a  certain  activity  may  be  triggered  only  when  the 
concentration  of  enabling  quantities  is  above  the  desired  threshold.  This  leads 
to  discrete  switching  between  active  and  dormant  states.  Second,  different 
models  may  be  appropriate  at  different  levels  of  concentration. 

Stochastic  behavior.  Evolution  of  entities  is  not  deterministic,  and  is  better 
captured  by  stochastic  models  that  allow  for  uncertainty  and  noise. 

These  characteristics  are  typical  of  high-level  models  of  embedded  software  such 
as  autonomous  communicating  mobile  robots.  For  describing  such  systems,  we 
have  developed  the  language  Charon  [2]  which  incorporates  ideas  from  con¬ 
currency  theory  (languages  such  as  CSP  [12]),  object-oriented  software  design 
notations  (such  as  Statecharts  [9]  and  UML  [3]),  and  formal  models  for  hybrid 
systems  (such  as  hybrid  automata  [1]  and  hybrid  I/O  automata  [15]).  The  key 
features  of  Charon  are: 

Architectural  hierarchy.  The  building  block  for  describing  the  system  ar¬ 
chitecture  is  an  agent  that  communicates  with  its  environment  via  shared 
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variables.  The  language  supports  the  operations  of  composition  of  agents  to 
model  concurrency,  hiding  of  variables  to  restrict  sharing  of  information,  and 
instantiation  of  agents  to  support  reuse. 

Behavior  hierarchy.  The  building  block  for  describing  flow  of  control  inside  an 
atomic  agent  is  a  mode.  A  mode  is  basically  a  hierarchical  state  machine,  that 
is,  a  mode  can  have  submodes  and  transitions  connecting  them.  Variables 
can  be  declared  locally  inside  any  mode  with  standard  scoping  rules  for 
visibility.  Modes  can  be  connected  to  each  other  only  via  well-defined  entry 
and  exit  points.  We  allow  sharing  of  modes  so  that  the  same  mode  definition 
can  be  instantiated  in  multiple  contexts.  Finally,  to  support  exceptions ,  the 
language  allows  group  transitions  from  default  exit  points  that  are  applicable 
to  all  enclosing  modes. 

Discrete  updates.  Discrete  updates  are  specified  by  guarded  actions  labeling 
transitions  connecting  the  modes.  Actions  can  have  calls  to  externally  defined 
Java  functions  which  can  be  used  to  write  complex  data  manipulations.  It 
also  allows  us  to  mimic  stochastic  aspects  through  randomization. 
Continuous  updates.  Some  of  the  variables  in  Charon  can  be  declared  ana¬ 
log ,  and  they  flow  continuously  during  continuous  updates  that  model  pas¬ 
sage  of  time.  The  evolution  of  analog  variables  can  be  constrained  in  three 
ways:  differential  constraints  (e.g.  by  equations  such  as  x  =  f(x,u,)),  alge¬ 
braic  constraints  (e.g.  by  equations  such  as  y  =  g(x,u)),  and  invariants  (e.g. 
\x  —  y\  <  e)  which  limit  the  allowed  durations  of  flows.  Such  constraints  can 
be  declared  at  different  levels  of  the  mode  hierarchy. 

Modular  features  of  Charon  allow  succinct  and  structured  description  of 
complex  systems.  Similar  features  are  supported  by  the  languages  Shift  [6]  and 
Stateflow  (see  www.mathworks.com).  In  Charon,  modularity  is  not  only  ap¬ 
parent  in  syntax,  but  we  are  developing  analysis  tools  (such  as  simulation)  that 
exploit  this  modularity.  Furthermore,  Charon  has  formal  foundations  support¬ 
ing  compositional  refinement  calculus  which  allows  relating  different  models  of 
the  system  in  mathematically  precise  manner.  A  formal  mathematical  descrip¬ 
tion  allows  us  to  develop  tools  for  computing  equilibria,  for  reachability  analysis 
and  for  analyzing  properties  like  stability  and  reachability. 

In  the  next  two  sections,  we  will  briefly  describe  case  studies  that  we  have 
used  to  investigate  the  hybrid  systems  approach  to  modeling  biological  systems, 
and  the  applications  of  Charon  to  these  systems.  We  will  also  illustrate  our 
approach  by  providing  preliminary  simulation  results. 


3  A  Repressilator  Network 

As  noted  in  [5],  most  biomolecular  systems  of  interest  involve  many  interactions 
connected  through  positive  and  negative  feedback  loops  and  an  understanding  of 
their  dynamics  is  hard  to  obtain.  In  this  section  we  will  describe  the  modeling  of 
a  specific  biomolecular  network.  We  will  model  a  repressilator  system  described 
in  [7].  First  we  provide  some  biological  background  information  and  describe  the 
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protein  network  used  in  [7],  and  then  describe  the  models  of  the  protein  network, 
including  examples  of  Charon  models.2 


3.1  The  Basic  Phenomena 

In  the  synthetic  oscillatory  network  described  in  [7],  networks  of  interacting 
biomolecules  carry  out  many  essential  functions  in  living  cells.  But  the  design 
principles  of  the  functioning  of  such  networks  still  remain  poorly  understood- 
even  in  relatively  simple  systems  [14].  The  authors  proposed  the  design  and 
construction  of  a  synthetic  protein  network  implementing  a  particular  function. 
Their  motivation  is  that  such  “rational  network  design”  may  lead  to  the  engi¬ 
neering  of  new  cellular  behaviors  and  to  improved  understanding  of  naturally 
occuring  networks. 

The  repressilator  system  described  in  [7]  contains  three  proteins,  namely 
lacl,  tetR,  and  cl.  The  protein  lacl  represses  the  protein  tetR,  tetR  represses  cl, 
whereas  cl  represses  lad,  thus  completing  a  feedback  system  called  a  repressi¬ 
lator  system.  The  dynamics  of  the  network  depend  on  the  transcription  rates, 
translation  rates,  and  decay  rates  of  proteins  and  messenger  RNAs.  Depending 
on  the  values  of  the  different  parameters  in  the  model,  the  system  might  converge 
to  a  stable  limit  cycle  or  become  unstable. 


3.2  Approaches  to  Modeling 

It  is  well  known  in  mechanics  and  thermodynamics  that  there  are  two  different 
approaches  to  modeling  systems  such  as  the  repressilator  system.  At  reasonably 
high  molecular  concentrations,  one  can  adopt  continuum  models  which  lend 
themselves  to  deterministic  models  involving  ordinary  and  partial  differential 
equations.  At  lower  concentrations,  the  discrete  molecular  interactions  become 
important  and  deterministic  models  are  difficult  to  obtain  [8]. 


The  Deterministic,  Continuous  Approximation.  We  will  consider  the 
three  repressor  protein  concentrations  pi,i  £  P  =  {lad,  tetR, cl}  and  their  cor¬ 
responding  mRNA  concentrations  rrii,i  £  P  as  continuous  dynamic  variables. 
The  system  kinetics  are  determined  by  the  following  six  coupled  first-order  dif¬ 
ferential  equations. 


drrii 

dt 

dpi 

dt 


—mi 


a 

l+p™ 


+  ao 


-(3(pi  -  mi) 


( i,j )  £  {(lacl, cl),  (tetR, lacl),  (cl, tetR)} 

2  For  more  information  on  Charon  or  sample  Charon  code,  please  check 
http://www.cis.upenn.edu/mobies/charon/  or  contact  ivancic@seas.upenn.edu. 
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The  equations  use  various  constants.  The  leakiness  of  the  promoter  a  is  the 
number  of  protein  copies  per  cell  produced  from  a  given  promoter  type  during 
continuous  growth  in  the  presence  of  saturating  repressor  amounts.  During  the 
absence  of  the  repressor,  we  have  a  +  «o  number  of  protein  copies  per  cell.  The 
ratio  of  the  protein  decay  rate  to  the  mRNA  decay  rate  is  denoted  by  /?,  while 
n  stands  for  the  so  called  Hill  coefficient. 

The  Stochastic,  Discrete  Approximation.  The  continuous  analysis  neglects 
the  discrete  nature  of  molecular  components  and  the  stochastic  character  of  their 
interaction  [7].  Following  [7],  we  adopt  the  stochastic  approximation  as  described 
by  Gillespie  in  [8].  The  various  proteins  and  mRNAs  are  modeled  by  discrete 
variables  corresponding  to  the  number  of  molecules  measuring  concentration, 
and  are  updated  at  discrete  time  intervals  by  stochastic  rules. 

3.3  Charon  Model 

In  this  section  we  will  present  the  repressilator  system  models  as  described  in  [7] 
using  the  Charon  language.  We  will  present  many  of  the  advantages  that  the 
Charon  language  has  to  offer  for  modeling  such  biomolecular  models. 

Our  model  will  define  a  generic  protein  model  as  an  agent  in  Charon.  We  will 
instantiate  this  agent  model  to  obtain  the  three  proteins  lad,  tetR,  and  cl.  The 
approximation  models  will  be  implemented  inside  the  modes  of  the  protein  agent. 
To  present  another  feature  of  our  language,  we  will  also  describe  a  combination 
of  the  discrete  and  the  continuous  model  into  one  modeling  system. 

The  Protein  Agent  in  the  Continuous  Approximation.  In  this  section 
we  will  describe  a  Charon  model  of  a  generic  protein  agent.  We  have  a  con¬ 
tinuous  input  variable  which  represents  the  repressor  protein  concentration  pr. 
This  means,  that  the  environment  of  this  protein  agent  supplies  the  value  of  this 
variable,  and  it  cannot  be  changed  by  the  protein  agent.  The  protein  agent  has  a 
continuous  private  variable  representing  the  messenger  R.NA  concentration.  Pri¬ 
vate  variables  cannot  be  seen  outside  the  agent  and  can  be  updated  internally  for 
internal  use  only.  The  output  of  the  protein  agent  is  a  continuous  variable  rep¬ 
resenting  the  protein  concentration.  Output  variables  are  updated  by  the  agent, 
and  can  be  used  as  input  variables  to  other  agents.  The  generic  protein  agent 
has  parameters  ao,  a,  j3,  n,po,  and  ?ti0-  By  instantiating  these  parameters  with 
values,  we  can  obtain  instantiated  protein  agents  representing  a  specific  protein. 
The  parameters  po  and  m o  will  be  used  for  initialization  purposes  and  stand  for 
the  initial  protein  concentration  and  the  initial  messenger  RNA  concentration 
respectively.  The  following  represents  the  corresponding  Charon  code. 

agent  contProtein  (real  pO  ,  mO  ,  alphaO  ,  alpha  ,  beta  ,  n){ 
write  analog  real  p  =  pO  ;  //protein  concentration 
read  analog  real  pR  ;  //repressor  protein  concentration 
private  analog  real  m  =  mO  ;  //messenger  RNA  concentration 
mode  cont  =  continuous  (  alphaO  ,  alpha  ,  beta  ,  n  )  ;  } 
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Fig.  1.  A  generic  protein  agent  for  the  continuous  approximation  model 


We  still  need  to  define  the  behavior  of  the  agent.  The  behavior  is  described 
by  the  modes  of  the  agent.  The  behavior  of  the  generic  protein  agent  is  defined 
in  cont,  which  is  an  instantiation  of  a  generic  continuous  mode  defined  by  the 
following  code.  A  graphical  version  of  the  generic  protein  model  can  be  found  in 
Figure  1. 

mode  continuous  (real  alphaO  ,  alpha  ,  beta  ,  n){ 
write  analog  real  p  ;  //protein  concentration 
read  analog  real  pR  ;  //repressor  protein  concentration 
private  analog  real  m  ;  //messenger  RNA  concentration 
diff  mRNA  {  d(m)  =  -m  +  alpha  /  (l+pR~n)  +  alphaO  } 
diff  proteinConcentration  {  d(p)  =  -beta  *  (p-m)  }  } 


Fig.  2.  Composed  repressilator  system  using  the  instantiated  generic  protein  agent 


Instantiation  and  Concurrency.  We  defined  a  generic  protein  agent  in  the 
previous  section.  We  have  to  instantiate  this  generic  agent  model  to  get  the 
three  proteins  used  in  the  system.  We  also  want  the  three  proteins  lad,  tetR, 
and  cl  to  run  in  parallel  and  to  influence  each  other.  Notice  the  use  of  renaming 
of  variables  to  couple  the  three  instantiated  protein  agents  to  influence  each 
other.  A  graphical  version  of  the  composed  system  is  illustrated  in  Figure  2.  The 
following  represents  the  corresponding  Charon  code  using  some  values  for  the 
parameters.  A  simulation  trace  generated  by  the  Charon  tool-set  is  given  in 
Figure  3. 
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agent  RepressilatorSystem  (){ 


private  analog  real  pi  , 

p2 

P3 

i 

agent  lacl  =  contProtein 

(  . 

.  ) 

[ 

P 

pR 

:=  pi 

p3  ] 

agent  tetR  =  contProtein 

(  . 

.  ) 

[ 

P 

pR 

:=  p2 

pi  ] 

agent  cl  =  contProtein 

(  . 

.  ) 

[ 

P 

pR 

:=  p3 

p2  ] 

} 

Fig.  3.  Simulation  trace  for  the  repressilator  system  showing  stable  oscillations  for  the 
three  protein  concentration  pi,p2,P3  over  time. 


The  Protein  Agent  in  the  Discrete  Approximation.  In  this  section  we 
will  present  a  possible  model  for  a  discrete  approximation  of  a  protein  agent.  As 
we  did  it  for  the  continuous  case,  we  will  again  define  a  generic  protein  agent, 
that  can  be  instantiated  to  build  a  system  of  proteins.  Our  model  works  as 
follows.  We  have  an  integer  variable  n  that  keeps  track  of  the  number  of  protein 
molecules  which  is  the  output  of  the  agent.  The  input  to  the  agent  is  the  number 
of  repressor  protein  molecules  hr.  Depending  on  various  parameters,  we  want 
to  increase  or  decrease  the  number  of  protein  molecules  by  one  at  a  time.  The 
basic  idea  is  to  use  stochastic  simulation  as  described  in  [8].  The  parameters 
that  influence  the  stochastic  simulation  are  binding  and  unbinding  of  proteins 
on  two-sided  promoters,  the  protein  and  mRNA  decay  rates,  and  translation. 
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protein  agent 
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discrete  mode 
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discreteTiming 

, 

continuous  mode 

*  . 

I« - 1 

Fig.  4.  A  generic  protein  agent  for  the  combined  framework  using  continuous  and 
discrete  approximation  model 


Combining  the  two  Models  into  one  Framework.  The  two  different  models 
for  the  repressilator  system  can  be  combined  into  one  framework.  The  basic  idea 
is  to  use  the  deterministic  continuous  model  whenever  the  concentration  of  the 
protein  is  high  enough,  whereas  we  would  switch  to  the  discrete,  stochastic  model 
if  the  concentration  would  fall  below  a  certain  threshold  value.  Figure  4  gives  an 
intuitive  graphical  representation  of  the  protein  agent  with  both  the  continuous 
and  discrete  approximation. 


4  Quorum  Sensing  in  Bacteria 

A  good  illustration  of  multicellular  behavior  in  prokaryotes  is  the  cell-density- 
dependent  gene  expression.  In  this  process,  a  single  cell  is  able  to  sense  when 
a  quorum  of  bacteria,  a  minimum  population  unit,  is  achieved.  Under  these 
conditions,  certain  behavior  is  efficiently  performed  by  the  quorum,  such  as  bio¬ 
luminescence,  which  is  the  best  known  model  for  understanding  the  mechanism 
of  cell-density-dependent  gene  expression.  In  this  section,  we  will  describe  a  hy¬ 
brid  system  model  that  captures  the  changes  in  dynamics  of  the  biochemical 
reactions  observed  in  the  literature  [13,16,17]. 


4.1  The  Basic  Phenomena 

Vibrio  fischeri  is  a  marine  bacterium  that  can  be  found  both  as  a  free-living 
organism  and  as  a  symbiont  of  some  marine  fish  and  squid.  As  a  free-living  or¬ 
ganism,  V.  fisheri  exists  at  low  densities  (less  than  500  cells  per  ml  of  seawater) 
and  appears  to  be  non- luminescent.  As  a  symbiont,  the  bacteria  live  at  high 
densities  and  are,  usually,  luminescent.  In  a  liquid  culture,  the  bacteria’s  level  of 
luminescence  is  low  until  the  culture  reaches  mid  to  late  exponential  phase.  A 
dramatic  increase  in  luminescence  is  observed  at  that  time  due  to  the  transcrip¬ 
tional  activation  of  the  lux  genes.  Once  the  bacteria  reach  stationary  phase,  the 
level  of  luminescence  decreases. 

The  /wx  regulon  [17]  contains  two  operons,  Ol  and  Or  (see  Figure  5).  The  left 
operon  Ol  contains  the  luxR  gene  encoding  the  protein  LuxR,  a  transcriptional 
activator  of  the  system.  The  right  operon  Or  contains  seven  genes  luxICDABEG. 
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Fig.  5.  A  portion  of  DNA  emphasizing  luxR  and  luxICDABEG  genes  and  the  binding 
sites  for  LuxR  complex  and  CRP 


Protein  LuxI,  the  product  of  the  luxl  gene  is  required  for  endogenous  production 
of  autoinducer ,  a  small  molecule  capable  of  diffusing  in  and  out  of  the  cell  mem¬ 
brane.  Genes  lux  A  and  luxB  encode  two  subunits  of  luciferase.  The  trio  luxC, 
luxD ,  and  luxE  code  for  the  subunits  of  a  protein  complex  which  provides  an 
aldehyde  substrate  for  luciferase.  The  function  of  luxG  is  unknown.  The  autoin¬ 
ducer  Ai  binds  to  protein  LuxR  to  form  a  complex  Co.  The  two  operons  are 
separated  by  a  regulatory  region  that  contains  a  binding  site  for  the  cyclic  AMP 
receptor  protein  CRP  and  a  binding  site  for  the  complex  Co. 

The  transcription  of  luxR  is  regulated  by  both  CRP  and  Co.  We  can  distin¬ 
guish  among  the  following  three  different  cases: 

—  Case  Oir  1  In  the  absence  of  the  autoinducer,  CRP  activates  Ol  expression 
by  initiating  two  R.NA  transcripts. 

—  Case  Ol~ 2  At  low  autoinducer  concentrations,  luxR  transcription  is  stimu¬ 
lated  by  increasing  CRP-dependent  transcription  and  by  Co-dependent  tran¬ 
scription  from  another  transcriptional  start  site. 

—  Case  Ol~ 3  At  high  autoinducer  concentrations,  luxR  transcription  is  re¬ 
pressed  through  a  second,  weaker  Co  binding  site  located  in  luxD. 

Likewise,  transcription  of  Or  is  regulated  by  both  CRP  and  Co.  We  distinguish 
two  different  cases: 

—  Case  Or-  1  In  the  absence  of  autoinducer,  CRP  represses  Or  transcription. 

—  Case  Or- 2  In  the  presence  of  autoinducer,  Co  activates  transcription  of  Or. 

These  cases  will  be  interpreted  as  modes  as  seen  later  in  the  paper. 

4.2  Mathematical  Model 

In  this  section,  we  develop  a  mathematical  model  for  the  luminescence  phe¬ 
nomenon  in  one  bacterium  of  V.  fischeri,  describing  the  concentrations  of  the 
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relevant  mRNA’s,  proteins,  and  small  molecules.  As  described  in  Section  4.1  the 
mechanism  of  transcription  activation  of  both  operons  is  highly  dependent  on 
the  concentration  of  autoinducer,  so  the  time  evolution  of  the  system  cannot  be 
described  by  one  set  of  continuous  differential  equations.  Combining  cases  for 
Ol  and  Or  given  in  the  previous  section,  yields  three  modes,  which  we  call  OFF, 
POS  and  NEG.  The  transitions  between  modes  are  governed  by  the  level  of  inter¬ 
nal  autoinducer  which  we  represent  by  [Ai] .  Mode  OFF  corresponds  to  very  low 

or  zero  concentration  of  autoinducer  ([Ai]  <  [Ai] _ )  within  the  bacterium  and  no 

luminescence  is  observed.  The  system  is  in  mode  POS  when  the  concentration 
of  internal  autoinducer  is  low  ( [ Ai]  _  <  [Ai]  <  [Ai]  , ) .  This  mode  corresponds 
to  positive  growth  and  increasing  concentration  of  autoinducer.  Luminescence 
is  observed,  as  are  higher  concentrations  of  proteins  LuxA,  LuxB,  LuxC,  LuxD, 
and  LuxE.  The  transition  to  mode  NEG  (negative  growth)  occurs  at  high  levels 
of  autoinducer  ([Ai]  >  [Ai]+). 

We  use  the  following  rate  equation  to  describe  the  concentration  for  any 
molecular  species  (mRNA,  protein,  protein  complex,  or  small  molecule)  [19]: 


d[x] 

——  =  synthesis  —  decay  ±  transformation  ±  transport 
at 


(i) 


The  synthesis  term  represents  transcription  for  mRNA  and  translation  for  pro¬ 
teins.  The  decay  term  represents  a  first  order  degradation  process.  The  transfor¬ 
mation  term  describes  reactions  such  as  cleavage  or  ligand-binding  that  do  not 
destroy  the  protein,  but  do  remove  its  ability  to  participate  in  specific  reactions. 
Finally,  molecular  species  may  participate  in  transport  processes,  like  passive 
diffusion  or  active  transport  through  a  membrane. 


The  biomolecular  system  can  be  described  in  a  nine  dimensional  state  space. 
The  nine  variables,  Xi,X2,  ■  ■  ■  ,xg,  describe  the  concentrations  of  different  mole¬ 
cules  as  follows: 


X\  =  mRNA  transcribed  from  Or, 

X2  =  mRNA  transcribed  from  Or, 

X3  =  protein  LuxR, 

X4  =  protein  LuxI, 

X5  =  protein  LuxA, 

Xq  =  protein  LuxB, 

X7  =  autoinducer  inside  the  bacterium  Ai, 

Xs  =  LuxR:Ai  complex  Co, 
xg  =  autoinducer  outside  the  bacterium  Aiex, 
where  Ai  is  the  dimensionless  version  of  [Ai]. 


For  simplicity,  we  have  assumed  that  the  concentrations  of  CRP  and  of  the 
substrate  necessary  for  endogenous  production  of  Ai  are  constant.  Further,  we 
have  neglected  the  decay  rates  for  chemical  compounds.  Finally,  we  assume  that 


3  In  [13],  the  differential  equations  for  the  low  autoinducer  concentration  are  described. 
The  model  presented  here  describes  a  wider  range  of  operating  conditions. 
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the  concentrations  of  LuxC,  LuxD,  and  LuxE  are  similar  to  those  of  LuxA  and 
LuxB. 

The  (continuous)  differential  equations  for  each  mode  are  of  the  form 
x  =  /*(x)  where  x  =  [*1,  x2,  ■  ■■,x9]T  €  IR9,  /*  =  [f{,  and 

i  G  {OFF,  POS,  NEG}.  The  components  of  the  vector  fields  are  explicitly 
given  by: 

f?FF  = 


nPOS 

J 1 


xNEG 

J 1 

fOFF 

J  2 

rPOS  _  rNEG 
J  2  —  J  2 

fi 

fi 

fi 

fi 

fi 


m 

4 


3c- 


V81 


si  ~4xl 


-J  I1X1 
-V2X2 


V2 


r,l/82 


-  X2 


82 


773  (xi  -  x3)  -  r37iAiX3X7 

774  (x2  —  X4)  —  r±X4 

775  (x2  -  X5) 

V6  (x2  -  x6) 

-777X7  +  r4x4  -  rmem  ( X7  -  Xg)  -  r3 7,rx3x7 


fi  =  -V8X8  +  r37'AiX3X7 

fi  =  -777X9  +  rmem{x7  -  Xg)  +  u 

where,  in  the  last  seven  equations  f!j  is  independent  of  the  mode.  All  the  quan¬ 
tities  in  the  above  model  are  non-dimensional,  rji  =  T3/F[i  where  To  is  the 
characteristic  time  constant  of  the  system  and  Hi  is  the  half-life  (inverse  of  the 
decay  rate)  of  molecule  x,.  vi:j  is  a  cooperativity  coefficient  while  Ky  describes 
the  potency  of  the  regulation  of  the  transcription  of  mRNA  j  by  protein  i.  r  de¬ 
notes  transformation  and  transfer  rates.  For  example  rmem  is  the  transfer  rate  of 
autoinducer  through  the  membrane  of  the  cell  while  7'37,r  and  r37^Ai  are  transfor¬ 
mation  rates  obtained  by  non-dimensionalizing  the  binding  rate  of  the  reaction 
between  Ai  and  LuxR  in  two  different  ways,  c  is  dependent  on  the  concentration 
of  CRP  and  its  affinity  to  the  corresponding  binding  site,  and,  as  stated  ear¬ 
lier,  is  assumed  to  be  constant.  Finally,  u  emulates  an  external  source  of  Ai  and 
is  used  to  simulate  the  sensitivity  of  the  bacterium  to  changes  of  autoinducer 
concentration  in  the  exterior. 

We  regard  u  as  an  input  to  our  system.  Since  proteins  LuxA  and  LuxB  are 
subunits  of  luciferase,  which  produces  luminescence,  it  is  reasonable  to  assume 
that  the  level  of  luminescence  is  proportional  to  the  product  of  the  concentrations 
of  LuxA  and  LuxB,  which  we  choose  to  be  the  output  of  the  system. 


4.3  Charon  Model 

The  behavioral  hierarchy  in  Charon  (see  Figure  6)  is  characterized  by  three 
different  behaviors  which  are  represented  by  three  different  modes,  namely  OFF, 


30 


R.  Alur  et  al. 


POS ,  and  NEG.  Many  of  the  differential  equations  governing  the  dynamics  of 
the  system  are  shared  between  the  modes.  We  will  introduce  the  notion  of  mode 
hierarchy  to  extract  the  shared  constraints.  Through  the  notion  of  submodes 
and  scoping,  we  can  simplify  the  description  of  the  respective  modes  OFF,  POS , 
and  NEG. 


Fig.  6.  Charon  structure  of  the  system 


Figure  7  illustrates  the  response  (i.e.,  luminescence)  of  the  bacterium  to  a 
perturbation  in  the  concentration  of  external  autoinducer  that  takes  the  form 
of  a  rectangular  pulse.  The  magnitude  of  the  step  has  been  chosen  to  make 
the  system  go  through  all  three  modes.  The  results  confirm  the  experimental 
observations  [17]:  luminescence  increases  during  mode  POS  and  decreases  in 
mode  NEG;  there  is  no  luminescence  in  mode  OFF.  The  switch  history  and  the 
time  evolution  of  the  concentrations  of  the  significant  molecules  in  the  system 
are  also  shown.  In  mode  OFF,  all  molecules  decay  to  zero,  except  for  mRNA 
Ol  and  the  corresponding  protein  R ,  as  expected.  For  a  short  time,  in  mode 
POS ,  all  the  concentrations  increase  until  the  internal  autoinducer  reaches  a  high 
concentration,  when  the  system  is  switched  to  mode  NEG.  In  this  last  mode, 
everything  decays  to  zero,  except  for  internal  autoinducer  which  can  reach  a 
stable  non-zero  value  dependent  on  the  size  of  the  step  of  external  autoinducer. 


5  Conclusions 

In  this  paper  we  have  shown  that  biological  cellular  networks  can  be  natu¬ 
rally  modeled  as  hybrid  systems.  In  particular,  the  protein  repressilator  system 
switches  between  a  continuous  deterministic  model  at  high  concentrations,  and 
a  timed,  discrete,  stochastic  model  at  low  concentrations.  Similarly,  the  lumi¬ 
nescence  control  of  Vibrio  fischeri  is  naturally  modeled  as  a  multi-modal  hybrid 
system,  resulting  in  simulations  that  are  in  accordance  with  experimental  obser¬ 
vations.  The  hybrid  nature  of  such  protein  networks  can  be  very  easily  expressed 
and  simulated  in  Charon,  which  may  offer  us  better  and  a  more  global  under¬ 
standing  of  biological  networks. 
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Fig.  7.  Increase  in  external  autoinducer  produces  luminescence:  (a)input  -  external 
source  of  autoinducer;  (b)  switch  history;  (c)  output  (luminescence)-  product  of  con¬ 
centrations  of  proteins  A  and  B-  (d)  and  (e)  time  -  evolution  of  concentrations; 


The  enormous  complexity  of  large  scale  biological  networks  will  present  us 
with  great  challenges  that  we  must  face.  Exploiting  the  structure  of  biological 
systems  will  be  critical  for  scaling  the  applicability  of  the  modeling,  analysis,  and 
simulation  tools.  It  is  therefore  extremely  encouraging  that  the  two  case  studies 
presented  in  this  paper  exhibit  the  architectural  paradigms  of  modern  software 
engineering. 

We  envision  the  link  between  hybrid  systems  technology,  and  biology  to 
strengthen.  The  scalable  nature  of  computational  tools  like  Charon  will  en¬ 
able  the  unified  and  improved  modeling  of  biological  cellular  networks,  leading 
to  better  understanding,  as  well  as  providing  us  with  the  opportunity  to  deter¬ 
mine  how  local  biological  changes  can  affect  global  behavior.  Conversely,  a  good 
understanding  of  the  robustness  of  noisy  biological  networks  will  lead  to  new 
approaches  to  designing  networked  embedded  systems. 

The  case  studies  also  highlight  the  need  for  developing  a  theory  of  stochastic 
hybrid  systems,  for  instance,  for  modeling  rate  equations  of  biochemical  reac¬ 
tions.  We  believe  that  mathematical  and  computational  tools  for  the  analysis 
of  such  systems  present  a  research  challenge  for  the  hybrid  systems  commu¬ 
nity,  while  presenting  a  significant  potential  for  greatly  impacting  post  genomics 
research. 
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